#### Aggregate Synth Diff-in-Diff ####

library(tidyverse)
library(ggpubr)
library(patchwork)
library(synthdid)
library(ggrepel)
library(xtable)
library(kableExtra)
library(knitr)
library(stargazer)
library(modelsummary)
library(purrr)
library(gsynth)
library(plotrix)
library(sandwich)
library(lmtest)

set.seed(7390)

trnt_agg <- read.csv('data/borough_trnt.csv')

### Including all Years and Councils ###

setup_agg <- panel.matrices(as.data.frame(trnt_agg),
                            unit = which(names(trnt_agg) == 'borough'),
                            time =  which(names(trnt_agg) == 'year'),
                            outcome = which(names(trnt_agg) == 'trnt'),
                            treatment = which(names(trnt_agg) == 'treat'))

tau.hat_agg <- synthdid_estimate(setup_agg$Y, 
                                 setup_agg$N0,
                                 setup_agg$T0)

tau.hat_agg

sdid_para_agg <- synthdid_plot(tau.hat_agg, facet.vertical=F,
                           control.name='Control', treated.name='Bromley',
                           lambda.comparable=TRUE, se.method = 'none',
                           trajectory.linetype = 1, line.width=2, 
                           effect.curvature=-.4, trajectory.alpha=.7, 
                           effect.alpha=1, diagram.alpha=1, 
                           onset.alpha=1) +
  theme(legend.position= 'bottom', 
        legend.direction='horizontal',
        legend.key=element_blank(), 
        legend.background=element_blank(),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        panel.grid = element_blank(),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        strip.text = element_text(size = 14))+
  scale_color_discrete('Borough')+
  scale_y_continuous(labels = scales::percent, 
                     limits = c(0.25,0.75), expand = c(0,0),
                     breaks = seq(0.3, 0.7, by = 0.1))+
  scale_x_continuous(breaks = seq(1966, 2018, by = 4))+
  coord_cartesian(xlim = c(1964,2018))+
  xlab('Year')+
  ylab('Turnout')

ggsave('plots/agg_para.png', sdid_para_agg, 
       width = 35, height = 25, units = 'cm')

### General Election Years Removed ###

setup_agg_gen <- panel.matrices(as.data.frame(trnt_agg %>% 
                                                filter(year != 1964, year != 1966,
                                                       year != 1970, year != 1974, year != 1979,
                                                       year != 1983, year != 1987, year != 1992,
                                                       year != 1997, year != 2001, year != 2005,
                                                       year != 2010)),
                                unit = which(names(trnt_agg) == 'borough'),
                                time =  which(names(trnt_agg) == 'year'),
                                outcome = which(names(trnt_agg) == 'trnt'),
                                treatment = which(names(trnt_agg) == 'treat'))

tau.hat_agg_gen <- synthdid_estimate(setup_agg_gen$Y, 
                                     setup_agg_gen$N0,
                                     setup_agg_gen$T0)
tau.hat_agg_gen

### Keeping only Tory Held Councils ###

trnt_agg_tories <- trnt_agg %>% 
  filter(borough == 'Barnet' | borough == 'Bromley' 
         | borough == 'Kensington and Chelsea' 
         | borough == 'Wandsworth' | borough == 'Westminster')


setup_agg_tories <- panel.matrices(as.data.frame(trnt_agg_tories),
                                   unit = which(names(trnt_agg_tories) == 'borough'),
                                   time =  which(names(trnt_agg_tories) == 'year'),
                                   outcome = which(names(trnt_agg_tories) == 'trnt'),
                                   treatment = which(names(trnt_agg_tories) == 'treat'))

tau.hat_agg_tories <- synthdid_estimate(setup_agg_tories$Y, 
                                        setup_agg_tories$N0,
                                        setup_agg_tories$T0)

### Creating Figure 6 ###

vc.hat.pl_agg <- vcov(tau.hat_agg, method='placebo')
se.hat.pl_agg <- sqrt(diag(vc.hat.pl_agg))
se.hat.pl_agg 

vc.hat.pl_agg_gen <- vcov(tau.hat_agg_gen, method='placebo')
se.hat.pl_agg_gen <- sqrt(diag(vc.hat.pl_agg_gen))
se.hat.pl_agg_gen

vc.hat.pl_agg_tories <- vcov(tau.hat_agg_tories, method='placebo')
se.hat.pl_agg_tories <- sqrt(diag(vc.hat.pl_agg_tories))
se.hat.pl_agg_tories

agg_mods <- ggplot()+
  geom_errorbar(aes(ymin =  
                      c(tau.hat_agg - 1.96*se.hat.pl_agg,
                        tau.hat_agg_gen - 1.96*se.hat.pl_agg_gen,
                        tau.hat_agg_tories - 1.96*se.hat.pl_agg_tories),
                    ymax = c(tau.hat_agg + 1.96*se.hat.pl_agg,
                             tau.hat_agg_gen + 1.96*se.hat.pl_agg_gen,
                             tau.hat_agg_tories + 1.96*se.hat.pl_agg_tories),
                    x = c('All Time-Periods',
                          'General Election Years Removed',
                          'Only Tory Held Councils')), 
                colour = '#F8766D',
                alpha = 0.8, width = 0.25, size = 1.5)+
  geom_point(aes(y = c(tau.hat_agg, tau.hat_agg_gen, 
                       tau.hat_agg_tories),
                 x = c('All Time-Periods',
                       'General Election Years Removed',
                       'Only Tory Held Councils')), 
             colour = '#00BFC4', alpha = 0.8, size = 2.5)+
  coord_cartesian(ylim = c(-0.15, 0.05))+
  scale_y_continuous(breaks = seq(-0.2, 0.1, by = 0.05),
                     labels = seq(-0.2, 0.1, by = 0.05), 
                     'Treatment Effect')+
  scale_x_discrete('Design')+
  theme_minimal()+
  theme(panel.grid = element_blank(), 
        axis.text = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 14),
        axis.title = element_text(size = 14))+
  geom_hline(yintercept = 0, colour = '#F8766D', linetype = 2)

agg_mods

ggsave('plots/agg_mods.png', agg_mods, 
       width = 35, height = 25, units = 'cm')


### Creating Table 3 ###

lambda_agg <- attr(tau.hat_agg, 'weights')$lambda
names(lambda_agg) <- colnames(setup_agg$W)[1:setup_agg$T0]
zeta.lambda_agg <- attr(tau.hat_agg, 'opts')$zeta.lambda_agg
time_weights_agg <- round(sort(lambda_agg, decreasing = T), 3)

time_weights_agg <- tibble(time_weights_agg) %>% 
  mutate(Year = names(time_weights_agg)) %>% 
  rename('Weight' = time_weights_agg) %>% 
  arrange(desc(Year)) %>% 
  relocate(Year)

time_weights_agg$Year <- as.character(time_weights_agg$Year)

print(xtable(time_weights_agg, type = 'latex', digits = 3, 
             label = 'tab:time_weights_agg', 
             caption = 'Time Weights '), include.rownames=FALSE,
      file = 'plots/time_weights_agg.tex')


### Creating Table 4 ###

omega <- attr(tau.hat_agg, 'weights')$omega
names(omega) <- rownames(setup_agg$W)[1:setup_agg$N0]
weights <- round(sort(omega, decreasing = T), 3)
zeta.omega <- attr(tau.hat_agg, 'opts')$zeta.omega

weights <- as.data.frame(weights)
weights <- cbind(Borough = rownames(weights), weights)
rownames(weights) <- 1:nrow(weights)

weights_agg <- knitr::kable(weights %>% 
                              rename('Weights'= weights), 
                            'pipe', digits = 3, 
                            caption = '\\label{tab:unit_weights_agg} Unit Weights for Borough Level Estimation')

save_kable(weights_agg, file = 'plots/weights_agg.Rmd')

